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ABSTRACT 



We show that numerical simulations of reionization that resolve the Lyman Limit systems 
(and, thus, correctly count absorptions of ionizing photons) have converged to about 10% level 
for 5 < z < 6.2 and are in reasonable agreement (within 10%) with the SDSS data in this 
redshift interval. The SDSS data thus constraint the redshift of overlap of cosmic Hii regions 
to zqvl — 6.1 ± 0.15. At higher redshifts, the simulations are far from convergence on the 
mean Gunn-Peterson optical depth, but achieve good convergence for the mean neutral hydrogen 
fraction. The simulations that fit the SDSS data, however, do not have nearly enough resolution 
to resolve the earliest episodes of star formation, and are very far from converging on the precise 
value of the optical depth to Thompson scattering - any value between 6 and 10% is possible, 
depending on the convergence rate of the simulations and the fractional contribution of PopIII 
stars. This is generally consistent with the third-year WMAP results, but much higher resolution 
simulation are required to come up with the sufficiently precise value for the Thompson optical 
depth that can be statistically compared with the WMAP data. 



Subject headings: cosmology: theory 
- galaxies: intergalactic medium 



cosmology: large-scale structure of universe - galaxies: formation 



1. Introduction 

Bad things comes in threes - but good things 
sometimes do too. At least, in twos. In the 
reionization research, the two are the latest anal- 
ysis of the Lyman-a absorption in the spectra of 
19 highest redshift SDSS quasars by Fan et al. 
(2006) and the large downward revision of the 
WMAP-mea.smed value for the Gunn-Peterson 
optical depth (Page 2006; Spergel 2006). 

The latest SDSS data are consistent with the 
end of reionization epoch at z ~ 6 (Fan et al. 
2006). The downward revision of the WMAP mea- 
surement now eliminates any need for complex sce- 
narios of early reionization (c.f. Melchiorri et al. 
2005, and references therein), which had to in- 
voke unknown or weakly constrained physics and 
which were at odds with the simplest reionization 
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scenario in which normal PopII stars are the dom- 
inant source of ionizing photons between z ^ 20 
and z ~ 5. 

In this simplest scenario (Gnedin 2000, 2004; 
Roy Choudhury & Ferrara 2006) - which we call 
the "Minimal Reionization Model" - PopIII stars 
and other exotic objects play an important but 
not a dominant role in reionizing the universe by 
z 6. The bulk of cosmic ionizations are pro- 
duced by normal PopII stars in sufficiently mas- 
sive galaxies (like the ones observed in UDF and 
GOODS). As these galaxies form, they create ion- 
ized (Hii) regions around them, which continue to 
expand. The universe is half-ionized (by volume) 
at z ^ 8 — 9, and shortly before z — 6 most of 
cosmic Hii regions merge during a relatively short 
period called "overlap" , completing the process of 
reionization of the universe. 

Since it now appears that the Minimal Reion- 
ization Model provides the best theoretical frame- 
work for the existing data, it would make sense 
to refine its predictions and expectations to make 
a more precise comparison with the SDSS and 
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Table 1: Simulation Parameters 



Run 










TT 


L4N128 


4 


1 


3.2 X lO'* 


6.2 


0.051 


L8N128 


8 


2 


2.6 X 10^ 


6.2 


0.048 


L8N256 


8 


0.64 


3.2 X 10'^ 


6.2 


0.056 



"Size of the computational box in Mpc. 
''Spatial resolution in comoving kpc. 
•^Mass resolution in Mq. 
''By construction. 



WMAP data. It is particularly timely now, since 

several different groups make simulations of reion- 
ization, and these simulations rarely agree with 
each other. 

In this paper wc revisit the agreement with the 
data (or lack of it) for the simulations reported in 
Gnedin (2004), as well as for the latest simulation 
that is substantially larger and has much higher 
spatial resolution. The existence of a set of simu- 
lations with varied box sizes and mass and spatial 
resolutions allows us not only to compare the data 
to the model, but to also estimate the level of nu- 
merical convergence of the simulations, and, thus, 
make a rough estimate of the current theoretical 
uncertainty. 

2. Simulations 

In this paper we use three simulations that have 
different spatial and mass resolutions and different 
sizes of the simulation volumes, which allows us to 
estimate the degree of numerical convergence for 
different physical quantities. All simulations have 
been performed with the Softened Lagrangian Hy- 
drodynamics (SLH) code (Gnedin 2000) using the 
Optically Thin Eddington Tensor (OTVET) ap- 
proximation (Gnedin & Abel 2001) for following 
the time-dependent and spatial-variable transfer 
of ionizing radiation in 3D. All three simulations 
use the same cosmological parameters as reported 
in Spergel et al. (2003) ^ 

The numerical parameters of the simulations 
are listed in Table 1. The first two simulations 
are ones labeled "A4" and "A8" in Gnedin (2004) 
and Fan et al. (2006), while the third one is a 
new simulation with 256^ dark matter particles, 




: 0.27, f2A = 0.73, h = 0.71, Ub = 0.04, ng = 1. 
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Fig. 1. A graphical representation of the spatial 
and mass resolutions for the three simulations we 
use in this paper. The horizontal axis is the spa- 
tial scale and the vertical axis is the mass scale. 
Spatial and mass scales resolved in the three sim- 
ulations are shown as three rectangles. 

the same number of baryonic quasi-Lagrangian 
cells, and about 1,300,000 stellar particles that 
have been forming continuously during the sim- 
ulation. In order to make references to specific 
simulations transparent, we label each simulation 
with a letter L followed by the value of the lin- 
ear size of the computational volume (measured 
in Mpc in comoving reference frame), followed 
by the letter N and the number of dark matter par- 
ticles (or baryonic cells) along one direction (128 
or 256). For example, L4N128 means a simula- 
tion with Ah~^ Mpc box size and with 128^ dark 
matter particles and an equal number of baryonic 
cells. 

Table 1 also gives the nominal spatial resolution 
of the simulations (as measured by the value of 
Plummer softening length - the real resolution be- 
ing a factor of 2-4 worse). Because of the specifics 
of the SLH method, a mesh with a larger number 
of cells can be deformed more, so the spatial res- 
olution of the large L8N256 simulation is, in fact, 
higher than the spatial resolution of the L4N128 
simulation, even if the mass resolutions of both 
simulations are identical. 

Finally, we also list in Table 1 values for the 
redshift of overlap, zqvl (Gnedin 2004), and the 
Thompson optical depth, tt, that we use below. 

It is convenient to represent the numerical res- 
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olution of a given simulation as a rectangle in the 
spatial scale vs mass scale plane, as shown in Fig- 
ure 1. The horizontal axis is the spatial scale {not 
a given spatial direction), and since a simulation 
has a finite spatial resolution and a finite box size, 
it is limited along the spatial axis at both ends. 
The same is true for the mass scale, so in the spa- 
tial scale - mass scale plane the simulation is rep- 
resented by a rectangle. The three simulations 
are shown in Fig. f . The L8N256 simulation has 
the largest spatial and mass dynamic range. The 
L8N128 run has the same box size as the L8N256, 
but lower spatial and mass resolution, while the 
L4N128 run has the same mass resolution as the 
L8N256 run, but lower box size (and somewhat 
lower spatial resolution, as explained above). 

Because of the computational expense, and be- 
cause our simulations become unreliable for 2; < 5 
(as we explain below), the L8N256 run has been 
continued only until z ^ 5. 

All simulations have been adjusted to best fit 
the mean transmitted flux data as explained in 
Gnedin (2004). 

3. How to Simulate Reionization Cor- 
rectly 

Before we can compare the simulations and the 
data, it is important to explain why our simula- 
tions are adequate for modeling reionization, de- 
spite the limited box size. Modeling reionization, 
after all, is all about counting absorptions of ion- 
izing photons correctly. It is well known that af- 
ter reionization, absorption of ionizing photons is 
dominated by the Lyman Limit systems (Miralda- 
Escude 2003). Obviously, the same should be true 
inside large enough Hii regions even during reion- 
ization, so resolving the Lyman Limit systems is 
crucial for counting the absorptions of ionized pho- 
tons correctly during and, perhaps, even before the 
overlap stage of reionization. 

To illustrate this point, we show in Figure 2 a 
sketch of Hii regions around two sources (a weak 
one and a strong one) before and after the overlap. 
Before overlap, most of photons emitted by a suf- 
ficiently weak source are absorbed by the neutral 
IGM just outside its I-front, while a sufficiently 
strong source, whose I-front reached the size com- 
parable to the mean free path of ionizing photons 
inside it, effectively reaches its "Stromgen sphere" , 



Before overlap After overlap 




Fig. 2. A sketch of two ionizing sources before 
and after overlap. Before overlap the mean free 
path for photons emitted by a weak source (in 
the upper left corner) is limited by the size of its 
Hii region (gray area), while the strong source (in 
the lower right corner) has reached its Stromgen 
sphere and all photons it emits are absorbed in 
Lyman Limit systems (small open cicles). After 
overlap, Lyman Limit system determine the mean 
free path for all sources. 

because inside its Hll region ionizations balance re- 
combinations^. After the overlap, all sources reach 
their "Stromgen spheres" , with ionization nearly 
balancing recombinations within the Lyman Limit 
systems (Miralda-Escude et al. 2000). 

The exact nature of Lyman Limit system is still 
rather poorly understood, but Lyman Limit sys- 
tems in a simulation can still be studied by closely 
mimicking the observational process. The simula- 
tions we use in this paper do indeed resolve the 
Lyman Limit system (Kohler & Gnedin 2006) - 
in fact, they do it only too well by overpredicting 
the observed numbers of Lyman Limit systems at 
2: ~ 4 by about a factor of 2 and by smaller factors 
at higher redshifts. The discrepancy is because 
the simulations do not properly include the ioniz- 
ing radiation from quasars, and they become in- 
adequate for describing the ionization state of the 
IGM for 2 < 5. This inadequacy at lower redshifts 
also becomes apparent in our results presented be- 
low. 

Based on the analysis of the Lyman Limit sys- 
tems from cosmological simulations we use here 
(Kohler & Gnedin 2006) it appears that spatial 

^Note, that we use the term "Stromgen sphere" here differ- 
ently from Shapiro & Giroux (1987), who did not account 
for Lyman Limit systems and thus concluded that cosmic 
Hii regions never reach their Stromgen spheres. 
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Fig. 3. — The moan Gunn-Pctcrson optical depth 
as a function of redshift for the simulations (gray 
lines) and data (black symbols). The dotted line 
shows the L4N128 run, the dashed line is for the 
L8N128 run, and the solid line is for our largest 
L8N256 run. In the latter case, the rms varia- 
tion around the mean is shown as a hatched band. 
Filled circles show the observational data from Fan 
et al. (2006) with both rms dispersion (thin error- 
bars) and errors of the mean (thick error-bars). 

resolution of at least 1-2 proper kpc is required to 
model the Lyman Limit systems correctly. Simu- 
lations of poorer spatial resolution would not be 
adequate for modeling cosmic reionization because 
they would underestimate absorption of ionizing 
radiation by a large factor. 

4. Results 

We first show in Figure 3 the evolution of the 
mean Gunn-Peterson optical depth, defined as 

TGP = -l0g((F(z))), (1) 

where {F[z)) is the mean (i.e. averaged over all 
possible directions) transmitted flux in the hydro- 
gen Lyman-a transition at a given redshift. It 
is important to underscore that the mean Gunn- 
Peterson optical depth is not an average of any- 
thing, and, following Fan ct al. (2006), we empha- 
size it by not using the overbar or the averaging 
operator () in its definition. 

The simulations are within about 10% of the 
data and within 10% of each other in the redshift 
interval 5 < z < 6.2, but the agreement with the 
data for the 5.8 < ^; < 6.2 interval is by construc- 



tion - we adjust the effective emissivity parameter 
in the simulations to fit the observed data, as is 
explained in Gnedin (2004). The agreement in the 
interval 5 < z < 5.8 is real - it reflects the fact that 
the ionization state of cosmic gas in the simula- 
tions is broadly consistent with the observational 
data (the following figures provide more detail on 
the level of agreement between the simulations and 
the data). 

The fact that simulations with various box sizes 
and resolutions differ from each by about 10% 
demonstrates that our numerical results have not 
converged to that level of precision. 

Formally, none of our three simulations is ac- 
tually consistent with the data at z < 5.8 in the 
statistical sense - the simulations lie within the rms 
fluctuation in the Gunn-Peterson optical depth for 
a range of redshifts, but outside the formal er- 
rors in the mean value of tgp, which are about 
5% for z < 5.8. However, the systematic errors 
in the measurement (continuum fitting and pos- 
sible contaminations from BAL/metal absorption 
in Lyman-a forest region) are likely to be of the 
order of 5 to 10%, so the agreement between the 
simulations and the data at 10% level is satisfac- 
tory at present and is an encouraging confirmation 
of the Minimal Reionization Model. 

We also emphasize the challenge facing the sim- 
ulations in the years to come. As the random 
and especially systematic errors in the data are re- 
duced, the current simulations will be hard pressed 
to fit the data at, say, 3 to 5% level. The future 
simulations will have to follow the radiative trans- 
fer correctly in fine detail to fit the observations, 
and this sensitivity to details offers a tremendous 
opportunities to learn about nature and distribu- 
tion of ionizing sources at z < 5.8. 

As we have mentioned above, the simulations 
become inadequate for z < 5 because they do 
not properly include the ionizing radiation from 
quasars. As the result, there are more Lyman 
Limit systems in the simulation (Kohler & Gnedin 
2006), the IGM in the simulations is more neutral, 
and the mean free path of ionizing radiation for 
z < 5 is shorter than arc actually observed. 

At z > 6.2 the simulations show a marked lack 
of convergence in the mean Gunn-Peterson optical 
depth. This difference cannot be explained away 
entirely by cosmic variance, since the rms fluctua- 



4 




Fig. 4. — The mean Gunn-Pctcrson optical depth 
as a function of redshift for the L8N256 simulation 
(gray lines) and data (black symbols). The solid 
line shows the Lyman-a measurement (the same 
line as in Fig. 3), and dashed and dotted lines 
show the Lyman-/? and Lyman-7 measurements 
respectively, rcscalcd to Lyman-a by factors 2.25 
and 4.4, as used by Fan et al. (2006). 

tions in, say, the L8N256 run are smaller than the 
difference between various runs. 

This difference, however, is not surprising. Af- 
ter all, a mean Gunn-Peterson optical depth of 
10 corresponds to the mean transmitted flux of 
only 1G~^ - in a simulation with 128"^ ^ 2 mil- 
lion cells only 20 transparent (F = 1) cells in the 
complete opaque {F = 0) medium would produce 
such a small mean transmitted flux. Thus, even 
discreetness effects in the simulation become im- 
portant at these low values of the mean transmit- 
ted flux, in addition to usual effects of limited res- 
olution, poorly known physics, lack of numerical 
convergence, etc - this is merely a statement of 
the fact that before the overlap any measurable 
transmission comes from the very tail of density 
and neutral fraction fluctuations, and, thus, is ex- 
ceedingly difficult to simulate correctly. The latter 
statement only applies to the transmitted flux in 
the hydrogen Lyman-a transition, other quantities 
like mean fractions or photoionization rates can be 
modeled more reliably before the overlap. 

It is also important to emphasize that the direct 
comparison between the simulations and the data 
for TGP > 10 is highly non-trivial. The observa- 
tional data constrain the Lyman-a optical depth 
directly only for tgp < 7, and, for higher val- 



Fig. 5. — The moan volume-weighted neutral hy- 
drogen fraction as a function of redshift for the 
simulations and data. The line and symbol mark- 
ings are as in Fig. 3. 

ues, properly rescaled constraints from Lyman-/3 
and Lyman-7 transitions are used. When we re- 
peat the Fan et al. (2006) rescaling procedure with 
our L8N256 simulation, we get a good agreement 
between the Lyman- /3 and Lyman-a for the scal- 
ing factor of 2.25, but get a higher optical depth 
for the Lyman-7 transition if we adopt Fan et al. 
(2006) scaling value of 4.4, while a smaller value 
of about 3.5 gives the best agreement between the 
Lyman-7 and Lyman-a measurements. 

Figure 5 shows the comparison between the 
SDSS data and the simulations for the mean vol- 
ume weighted neutral fraction^ function of 
redshift. In this quantity, the convergence of the 
simulations is much better, and the data are con- 
sistent (at 10 to 20% level) with the simulation 
results. This agreement is not surprising, given 
the agreement in the mean Gunn-Peterson opti- 
cal depth - the neutral fraction (and the mean 
free path that we discuss below) is not measured 
directly from the observations, but is rather de- 
rived from the Gunn-Petrson optical depth mea- 
surement based on assumed density distribution 
of the IGM and photoionization equilibrium. 

In Figure 6 we show the evolution of the mean 

*The comparison is much more difficult for tlic mean mass 
weighted neutral fraction, because the simulations do not 
resolve Damped Lyman-a systems, that are known to con- 
tain most of neutral gas in the universe at z < 4, and are 
likely to contain at least a substantial fraction of all neutral 
gas at 2 < 6 as well. 
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Fig. 6. The mean free path of ionizing radia- 
tion as a function of redshift for the simulations 
and data. The hne and symbol markings are as in 
Fig. 3. The filled black square at z = 4 shows the 
observational determination of the mean free path 
by Miralda-Escudc (2003), and the hatched black 
band is the extrapolation of the observed evolu- 
tion of the Lyman Limit systems from Storrie- 
Lombardi et al. (1994) to ^; > 4. 

free path of the Lyman Limit photons in the sim- 
ulations and in the SDSS data, as well as a plau- 
sible extrapolation to high redshift of the mean 
free path from the Lyman Limit systems alone 
(Storrie-Lombardi et al. 1994; Miralda-Escude 
2003). As we have mentioned before, the simu- 
lations slightly underestimate the mean free path 
because they overpredict the abundance of the 
Lyman Limit systems up to a factor of 2 (Kohler 
& Gnedin 2006). The SDSS data only go up to 
the overlap epoch at z « 6.2, and are not yet able 
to probe the approach to the overlap at z > 6.2. 

The SDSS estimate of the mean free path lies 
somewhat below the extrapolation from the Ly- 
man Limit system measurement at z = 4. This 
difference is most likely not significant, as the ex- 
trapolation assumes a power-law evolution of Ly- 
man Limit systems in redshift, which is too simple 
to be precise. But the general agreement between 
the extrapolation and the SDSS result indicates 
that the mean free path of Lyman Limit photons 
is limited by the Lyman Limit systems for z < 5.8, 
as is generally believed. 

In the simulations, the strong deviation from 
the power-law evolution of the mean transmitted 
path with redshift at ^; > 6 indicates the over- 



Fig. 7. The Thompson optical depth as a func- 
tion of resolution for the L8N128 and L8N256 
simulations (black squares). The gray and black 
hatched bands show the constraints on tt from 
the WMAP polarization measurement alone (tt = 
0.10 ± 0.03, Page (2006)) and from the combined 
WMAP+SmS data (tt = 0.08 ± 0.03, Spergel 
(2006)) respectively. Black dotted, dashed, and 
solid curved lines show the fits to the simulation 
results in the form T'y{N) — Tr^ + A/N*^ for a = 1, 
1/2, and 1/4 respectively, while horizontal lines of 
the same type give the converged values Too- 

lap stage of reionization - during which the mean 
free path is still determined by a typical size of 
Hll regions rather than by the Lyman Limit sys- 
tems. If we take the agreement between the SDSS 
data and the simulations for z < 6.2 as an indica- 
tion that the rapid decrease in the SDSS estimate 
of the mean transmitted path indeed corresponds 
to the overlap of cosmic Hll regions, then we can 
adopt the last SDSS data point as a plausible es- 
timate of the time of overlap (formally defined as 
the moment when the mean free path grows most 
rapidly), giving zqvl = 6.1 ± 0.15. 

In particular, the SDSS data limit the typical 
size of a cosmic Hil region to less than about lh~^ 
proper Mpc at z 6.2. 

Finally, we can compare the optical depth to 
Thompson scattering, tt, to the WMAP measure- 
ments. These values for the three simulations are 
given in Table 1 and also shown in Figure 7. As 
one can see, these values arc far from being con- 
verged, because the smallest objects that form first 
(and are missed in finite resolution simulations) 
can contribute significantly to the Thompson op- 



6 



tical depth. In numerical analysis it is customary 
to estimate the converged result for any physical 
quantity C from a set of simulations with var- 
ied mesh size N by fitting the computed values 
of C{N) with a functional form 

C{N)=C^ + ^. (2) 

In order to determine three coefficients Cqo, A, 
and a, three simulations with three different mesh 
sizes are required. Unfortunately, in our case it 
is not possible to get another simulation with a 
significantly different value of N: running a 512^ 
SLH simulation is completely unfeasible now, and 
a 64^ simulation is so small that it does not fit 
the SDSS data at all. Also, using a value of N 
somewhere between 128 and 256 would not help, 
as Fig. 7 illustrates. 

In order to illustrate the inadequacy of our sim- 
ulations for giving an accurate estimate of tt, wc 
fitted equation (2) for the L8N128 and L8N256 
simulations (since the resolution is the most im- 
portant quantity for modeling the highest rcdshift 
star formation) keeping a as a free parameter, 
and fits for a = 1, 1/2, and 1/4 are shown in 
Fig. 7. As anyone can sec, the simulations arc 
non-conclusive; depending on the rate of numeri- 
cal convergence, the simulations can produce any 
value for between about 6 and 10%. 

In addition, a significant (but not dominant) 
contribution from PopIII stars is also possible. 

5. Conclusions 

We have shown that numerical simulations of 
reionization that resolve the Lyman Limit systems 
(and, thus, correctly count absorptions of ionizing 
photons) are close (within 10% or so) to the SDSS 
data for a variety of measured and deduced quan- 
tities for z < 6.2. The SDSS data thus constraint 
the redshift of overlap to zqvl = 6.1 ± 0.15. 

The mean free path of ionizing photons is lim- 
ited by the Lyman Limit systems at ^; < 6.0, and is 
smaller at higher redshifts, reflecting a finite size of 
a typical cosmic Hil region (less than lh~^ proper 
Mpc at ^ = 6.2). 

Our simulations, however, do not have nearly 
enough resolution to resolve the earliest episodes 
of star formation, and are very far from converg- 
ing on the precise value of the optical depth to 



Thompson scattering - any value between 6 and 
10% is possible, depending on the convergence rate 
of the simulations and the fractional contribution 
of PopIII stars. This is generally consistent with 
the third-year WMAP results, but much higher 
resolution simulation are required to come up with 
the sufficiently precise value for the Thompson op- 
tical depth that can be statistically compared with 
the WMAP data. 

While our simulations agree with the SDSS 
data within about 10% for 5 < 2: < 6.2, the level of 
numerical convergence of simulations in this red- 
shift interval is 10% at best. As the data improve 
(mostly by reducing the systematic errors), the 
constraining power of the SDSS data will be suffi- 
cient to place non-trivial demands on the degree of 
realism of any simulation that attempts to fit the 
data statistically. Not only the cosmological pa- 
rameters in the simulation must be precise enough, 
but it is likely that fine details of radiative trans- 
fer (the relative roles of quasars and galaxies, their 
spatial clustering, accuracy of numerical schemes, 
etc) have to be done accurately as well. These 
requirements will present a challenge and a moti- 
vation for the future theoretical work on modeling 
cosmic reionization with greater precision than has 
been possible so far. 
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